{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "dd6bb138",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "e8f9f7ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Program_3():\n",
    "    ############## Parameters\n",
    "    global c, xmu, eps0, asize\n",
    "    global nx, ny, nt, mxst, mxnd, myst, mynd\n",
    "    global dt, ds\n",
    "    global Ez, Hx, Hy #Create E and H field components.\n",
    "    # global data type enables the global scope of the variables\n",
    "    # within the code.Unlike global variables they can not be accessed\n",
    "    # outside the code\n",
    "    global mediaEz, mediaHx, mediaHy\n",
    "    global Ca, Cb, Da, Db # Define material based coefficients\n",
    "    global c1,c2,c3,c4,c5\n",
    "    global Ez1, Ez2, Ez3, Ez4, Ez5 # For Bubbling of E-fields in Liao ABC\n",
    "    global ABC_order # Order of the Liao ABC \n",
    "    global strip\n",
    "    \n",
    "    c = 2.99792458e8;\n",
    "    xmu = 4*math.pi*1e-7;\n",
    "    eps0 = 8.854187817e-12;\n",
    "    asize = 5; # Space Dimension in meters\n",
    "    nx = 80;      # Number of cells in x-direction\n",
    "    ny = 100;     # Number of cells in y-direction\n",
    "    nt = 400;     # Number of time steps\n",
    "    mxst = 17;    # Start of PEC section in x-direction\n",
    "    mxnd = 49;    # End of PEC section in x-direction\n",
    "    myst = 33;    # Start of PEC section in y-direction\n",
    "    mynd = 65;    # End of PEC section in y-direction\n",
    "    strip=30\n",
    "    \n",
    "    Ez = np.zeros((nx, ny)) # z-component of E-field\n",
    "    Hx = np.zeros((nx, ny)) # x-component of H-field\n",
    "    Hy = np.zeros((nx, ny)) # y-component of H-field\n",
    "\n",
    "    mediaEz = np.ones((nx, ny)) # z-component of E-field\n",
    "    mediaHx = np.ones((nx, ny)) # x-component of H-field\n",
    "    mediaHy = np.ones((nx, ny)) # x-component of H-field\n",
    "\n",
    "    Ca = np.zeros((2, 1)) # x-component of H-field\n",
    "    Cb = np.zeros((2, 1)) # x-component of H-field\n",
    "    Da = np.zeros((2, 1)) # x-component of H-field\n",
    "    Db = np.zeros((2, 1)) # x-component of H-field\n",
    "\n",
    "    Ez1 = np.zeros((nx,ny));\n",
    "    Ez2 = np.zeros((nx,ny));\n",
    "    Ez3 = np.zeros((nx,ny));\n",
    "    Ez4 = np.zeros((nx,ny));\n",
    "    Ez5 = np.zeros((nx,ny));\n",
    "    \n",
    "    ds = asize/(mxnd - mxst - 1); # Length Increment\n",
    "    dt = ds/(c*math.sqrt(2)); # Stability Condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "85db743d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def define_media(iflaga):\n",
    "    global nx, ny, mxst, mxnd, myst, mynd\n",
    "    global mediaEz, mediaHx, mediaHy\n",
    "    if (iflaga == 2):\n",
    "        for i in range(nx):\n",
    "            for j in range(ny):\n",
    "                if (i >= mxst-1 and i <= mxnd-1):\n",
    "                    if (j >= myst-1 and j <= mynd-1):\n",
    "                        mediaEz[i, j] = 2\n",
    "        for i in range(nx):\n",
    "            for j in range(ny):\n",
    "                if (i >= mxst-1 and i <= mxnd-1):\n",
    "                    if (j >= myst-1 and j <= mynd - 2):\n",
    "                        mediaHx[i, j] = 2\n",
    "        for i in range(nx):\n",
    "            for j in range(ny):\n",
    "                if (i >= mxst and i <= mxnd - 1):\n",
    "                    if (j >= myst-1 and j <= mynd-1):\n",
    "                        mediaHy[i, j] = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "db398631",
   "metadata": {},
   "outputs": [],
   "source": [
    "def define_coefficients():\n",
    "    global Ca, Cb, Da, Db\n",
    "    global xmu, eps0, dt, ds\n",
    "    dte = dt/(ds*eps0)\n",
    "    dtm = dt/(ds*xmu)\n",
    "    Da[0] = 1\n",
    "    Db[0] = dtm\n",
    "    Ca[0] = 1\n",
    "    Cb[0] = dte\n",
    "    Da[1] = 0\n",
    "    Db[1] = 0\n",
    "    Ca[1] = 0\n",
    "    Cb[1] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "78c0bc5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def define_Liao():\n",
    "    global c1,c2, c3, c4, c5 ; # Define material based coefficients\n",
    "    global ABC_order # Order of the Liao ABC \n",
    "    if ABC_order == 5:\n",
    "        c1=5;\n",
    "        c2=10;\n",
    "        c3=10;\n",
    "        c4=5;\n",
    "        c5=1;\n",
    "    elif ABC_order == 4:\n",
    "        c1 = 4;\n",
    "        c2 = 6;\n",
    "        c3 = 4;\n",
    "        c4 = 1;\n",
    "        c5 = 0;\n",
    "    elif ABC_order == 3:\n",
    "        c1 = 3;\n",
    "        c2 = 3;\n",
    "        c3 = 1;\n",
    "        c4 = 0;\n",
    "        c5 = 0;\n",
    "    else:\n",
    "        print('Error: Wrong Value')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "ed7d3515",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Liao_ABC():\n",
    "    global c1, c2, c3, c4, c5 ; # Define material based coefficients\n",
    "    global Ez ; # E and H field components.\n",
    "    global Ez1, Ez2, Ez3, Ez4, Ez5 # For Bubbling of E-fields in Liao ABC\n",
    "    global nx,  ny  \n",
    "    global ABC_order\n",
    "    for j in range(ny):\n",
    "        Ez[1,j]=c1*Ez1[2,j]-c2**Ez2[3,j]+c3*Ez3[4,j]-c4*Ez4[5,j]+c5*Ez5[6,j] # left side\n",
    "    for j in range(ny):\n",
    "        Ez[nx-1,j]=c1*Ez1[nx-2,j]-c2*Ez2[nx-3,j]+c3*Ez3[nx-4,j]-c4*Ez4[nx-5,j]+c5*Ez5[nx-6,j]\n",
    "    for i in range(2,nx-1):\n",
    "        Ez[i,1]=c1*Ez1[i,2]-c2*Ez2[i,3]+c3*Ez3[i,4]-c4*Ez4[i,5]+c5*Ez5[i,6]\n",
    "    for i in range(2,nx-1):\n",
    "        Ez[i,ny-1]=c1*Ez1[i,ny-2]-c2*Ez2[i, ny-3]+c3*Ez3[i,ny-4]-c4*Ez4[i,ny-5]+c5*Ez5[i,ny-6] #top\n",
    "    \n",
    "    if ABC_order == 5:\n",
    "        Ez5 = Ez4;\n",
    "        Ez4 = Ez3;\n",
    "        Ez3 = Ez2;\n",
    "        Ez2 = Ez1;\n",
    "        Ez1 = Ez;\n",
    "    elif ABC_order == 4:\n",
    "        Ez4 = Ez3;\n",
    "        Ez3 = Ez2;\n",
    "        Ez2 = Ez1;\n",
    "        Ez1 = Ez;\n",
    "    elif ABC_order == 3:\n",
    "        Ez3=Ez2;\n",
    "        Ez2=Ez1;\n",
    "        Ez1=Ez;\n",
    "    else:\n",
    "        print('Error: Wrong Value')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "7b31e456",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Source(n, sources):\n",
    "    global Ezs\n",
    "    if sources == 2:\n",
    "        xndec = 10.0\n",
    "        xn0 = 4*xndec\n",
    "        Ezs = math.exp(-((n-xn0)/(xndec))**2)\n",
    "    elif sources == 1:\n",
    "        if ( n >=1 and n <= 10):\n",
    "            Ezs = math.sin((n)*math.pi/10)\n",
    "    return Ezs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "d3629e3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def adv_Ez(n, sources):\n",
    "    global Ez, Hx, Hy\n",
    "    global mediaEz\n",
    "    global Ca, Cb\n",
    "    global nx, ny\n",
    "    for i in range(nx):\n",
    "        for j in range(ny):\n",
    "            m  = int(mediaEz[i, j]-1)\n",
    "            if (i == 5):\n",
    "                Es= Source(n, sources)\n",
    "            else:\n",
    "                Es= 0     \n",
    "            if (i >= 1 and j >=1):\n",
    "                Ez[i, j] = Ez[i, j]*Ca[m] + Cb[m]*(Hy[i, j] - Hy[i-1, j] - (Hx[i, j] - Hx[i, j-1]))+Es\n",
    "            elif (i>=1 and j == 0):\n",
    "                Ez[i, j] = Ez[i, j]*Ca[m] + Cb[m]*(Hy[i, j] - Hy[i-1, j]- Hx[i, j])+Es"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "5fbbd083",
   "metadata": {},
   "outputs": [],
   "source": [
    "def adv_H(n):\n",
    "    global Ez, Hx, Hy\n",
    "    global mediaHx, mediaHy\n",
    "    global Da, Db\n",
    "    global nx, ny\n",
    "    for i in range(nx):\n",
    "        for j in range(ny-1):\n",
    "            m = int(mediaHx[i, j]-1)\n",
    "            Hx[i, j] = Hx[i, j]*Da[m] - Db[m]*(Ez[i, j+1] - Ez[i, j])\n",
    "    for i in range(nx-1):\n",
    "        for j in range(ny):\n",
    "            m = int(mediaHy[i, j]-1)\n",
    "            Hy[i, j] = Hy[i, j]*Da[m] + Db[m]*(Ez[i+1, j] - Ez[i, j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "0092909d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_surface_plot(field):\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "    xd = np.linspace(0, 10, 100)\n",
    "    yd = np.linspace(0, 10, 80)\n",
    "    [xdg, ydg] = np.meshgrid(xd, yd)\n",
    "    dem3d = ax.plot_surface(xdg, ydg, field, cmap='rainbow', edgecolor='none')\n",
    "    fig.colorbar(dem3d, shrink=0.5, aspect=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "fbfa274b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_line_plot(Ez_5, Ez_35, Ez_65, Ez_95,iflaga, source):\n",
    "    fig = plt.figure()\n",
    "    plt.subplot(411)\n",
    "    plt.plot(range(len(Ez_5)),Ez_5,'k')\n",
    "    plt.xticks([])\n",
    "    plt.legend(title='n=5', loc='best', frameon=False)\n",
    "    plt.subplot(412)\n",
    "    plt.plot(range(len(Ez_35)),Ez_35, 'k')\n",
    "    plt.xticks([])\n",
    "    plt.legend(title='n=35',loc='best', frameon=False)\n",
    "    plt.subplot(413)\n",
    "    plt.plot(range(len(Ez_65)),Ez_65,'k')\n",
    "    plt.xticks([])\n",
    "    plt.legend(title='n=65',loc='best', frameon=False)\n",
    "    plt.subplot(414)\n",
    "    plt.plot(range(len(Ez_95)),Ez_95,'k')\n",
    "    plt.xticks([])\n",
    "    plt.legend(title='n=95',loc='best', frameon=False)\n",
    "    if iflaga == 1 and source == 1:\n",
    "        fig.suptitle('Line Plots for E-field with no obstacle for a sinusoid source')\n",
    "    elif iflaga == 2 and source == 1:\n",
    "        fig.suptitle('Line Plots for E-field with a PEC box for a sinusoid source')\n",
    "    elif source == 2:\n",
    "        fig.suptitle('Line Plots for E-field with no obstacle for a Gaussian source')\n",
    "    # plt.title('Line Plots for E-field with no obstacle for a sinusoid source')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "ca72947c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def zeroing():\n",
    "    global nx, ny\n",
    "    global Ez, Hx, Hy\n",
    "    global mediaEz, mediaHx, mediaHy\n",
    "    global Ca, Cb, Da, Db\n",
    "\n",
    "    Ez = np.zeros((nx, ny))\n",
    "    Hx = np.zeros((nx, ny))\n",
    "    Hy = np.zeros((nx, ny))\n",
    "    \n",
    "    Ez1 = np.zeros(nx,ny);\n",
    "    Ez2 = np.zeros(nx,ny);\n",
    "    Ez3 = np.zeros(nx,ny);\n",
    "    Ez4 = np.zeros(nx,ny);\n",
    "    Ez5 = np.zeros(nx,ny);\n",
    "\n",
    "\n",
    "    mediaEz = np.ones((nx, ny))\n",
    "    mediaHx = np.ones((nx, ny))\n",
    "    mediaHy = np.ones((nx, ny))\n",
    "\n",
    "    Ca = np.zeros((2, 1))\n",
    "    Cb = np.zeros((2, 1))\n",
    "    Da = np.zeros((2, 1))\n",
    "    Db = np.zeros((2, 1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "192f4894",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No handles with labels found to put in legend.\n",
      "No handles with labels found to put in legend.\n",
      "No handles with labels found to put in legend.\n",
      "No handles with labels found to put in legend.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsc0lEQVR4nO3de3xU9Z3/8dcnCQlX5RIyQggEbywQJUqUrTfogqDUeilt11pv0B9qK+tulbTaYtdWUSvWWrWW6mq7y6bSqlVbpN7axctahaSyCCpgJZlAhAQwAhIjSb6/P85JGCYzuc9kmHk/H495ZOac7znf75lzct5z7uacQ0REUldabzdARER6l4JARCTFKQhERFKcgkBEJMUpCEREUpyCQEQkxfVKEJjZmWa2sRfqLTezGXGo55tmtsPM9pnZsFjXF6H+081ss1//hWb2JzO7ooPDOjM7Nkq/K83stW62rc15b2b5fhsyulNPTzOzX5vZbb3chi4vv729TCYaM/u6mb3Q2+1IFDENgmgLrnPuVefcuBjV6czsE3+B32Zm95hZeifHMc3Mtnax/j7APcBM59xA59yurownbJzlZlbnT1Pz64E2BvkR8IBf/9POuXOdc//Z3Xb0hPB5H69w7i0JEiA9vkx2oQ3HmdlyM6sxsz3+D5X7zWxUvNsC4Jwrcc7N7I26E1Gy7hqa5JwbCEwHLgHmx7HuANAX2NDZAc0TbZ580f8nbn4taGNUY7pSvyStWC2THR3HscCbQBVwknPuCOB04O/AGd0ZdyrqiXnSinMuZi+gHJgRofs0YGtYuYXAOuBj4LdA35D+5wFrgVrgdeDENup0wLEhnx/H+3V8SHuALOBevIWzyn+fBQwA6oAmYJ//GgmcCpQCe4AdwD0R6j4e+MRvwz7gL37304A1/rStAU4LGWYVsBj4X7/eYzv6PUaZ/r/7ba/z25Dl1/H/QsrMA94FPgKeB8ZE+v6AYcAf/GleDdwKvBal3v8EbvDf5/rj+Zb/+VhgN2Ch8x5YFtbW7wD5/rBXAEFgJ/D9Nqb318DPgWeBvXgrnGNC+kf97iOMa7z/XdXirTTPD6tnKfCiX8/Lzd+bP10/Bar9etYBBcBVwAHgM3/6/uiXv9GfT3uBd4CLwtox358/zf1PjrD8poWMZxfwO2BoHJfJNqchrOx/N097G2WGACuAGrzlcgUwKtr/AHAL8N/++75+Hbv8ebcGCPj9rgQ+8Nu5Bfh6SPfXQsb3M6ASb1kvA84Mq+t3wH/549kAFEWZjojLgt/vSH8cNUAFsAhIC58e/3O+P88yos0TYCLe8rgbb530vc4sG4e0uyMrl66+wmdeSPdptA6C1Xgr3KF4/wTX+P1O9r/UKUA63gqiHMiKUmfoimwCsB34RoR/pB8BbwA5wHC8gLk1Uvv8bn8FLvPfDwT+MUr94TNwKN6CfRmQAXzN/zwsZAYH/ZmaAfTp6PfY0e+dkCAALgTex1vpZfgL4+tRvr/l/kI0AG/Fto3oQTCPgyu6S/yF8Lch/Z5pY96HtrX5+3sY6AdMAuqB8VHq/TXeP8Kp/vSUAMs78t2HjaeP/718D8gE/gnvn35cSD17gbPwwvVnzd8FMAtv5TEYb0UwHhgRMtxtYXV9BW9ZTwP+GW9FPSKk3zbgFH9cx3IwcFq+K+Df8JbfUX57fgk8FsdlMuo0RCi7HbiynWV2GDAH6A8MwvsB93Qby8ktHAyCq4E/+sOmA5OBI/CW2z0h83AEMNF/fyWHBsGlfhsygBv8NvcNqetTYLY//juAN6JMR1vLwn8Bz/jTlw9s4uC6qWV6osyz8HkyCPjQb2tf//OUzi4bLfV1dOXSlVf4zAvpPo3WK4NLQz7fBSz13/8CfwUd0n8jMDVKnc6f+R/hrYxu42DqtrTH7zc7bAaWR2qf3+0V4IdAdjvTHD4DLwNWh5X5K/4/hj+Df9SB73Ef3q+d5tf8jn7vHBoEf2pe+PzPacB+Dq5sHN7KJx3v1+w/hJS9nehBcIzfrjS8X85Xc/CX/38C17cx7yMFQeivwdXAxVHq/TXwHyGfZwPvdeS7D+t+Jt4/f1pIt8eAW0LqWR7SbyDQCOThhcYm4B9Dhw8Z7rZIbQ8psxa4wH//PPCv7c1XvB9L00P6jfDnV0Y8lsm2piFCvwbgnJDPC/xlZR/wcJRhCoGP2lhObuFgEMwjwp4CvCCoxQuYfmH9riTKsuz3/whvF3NzXS+F9JsA1EUZLuKygPf/VA9MCOl2NbAqfHqizLND5gleeL8VpQ0dXjaaX4l0jGB7yPv9eP9o4O3vvsHMaptfeP98I9sY18nOuSHOuWOcc4ucc00RyozE2zxrVtHOOL+Bt5n9npmtMbPz2pmeaPU015Ub8rmyA+O50Dk3OOT1MICZbQg5gHxmB8YzBvhZyHfZvMsmN6zccLxfHqFtC5+OFs65v+P9YxfirVRXAFVmNg6YircrpTOiLQ+dKduR756QspVhy0rU+eSc24f33Y10zv0FeABvF9UOM3vIzI6I1lgzu9zM1obMgwIg2++dh/cjpT1jgKdCxvEuXjAFOjBst5fJdqYh3C68lREAzrkHnHOD8XbH9vHH19/MfmlmFWa2B++H1+AOnuixDC9Al5tZlZndZWZ9nHOf4G2tXAN8aGbPmtk/RJmeG8zsXTP72J+eI8OmJ3wZ6xvpzLY2loVsvC3N8HVOpGUxmtB50tZy0ullI5GCIJpKYHHYSrC/c+6xbo63Cu8Lazba7wZeEh/CObfZOfc1vF1JPwaeMLMBXainua5toaPvaKMjtGuiO3gA+dUODFIJXB32ffZzzr0eVq4G75dcXli72/Iy8GUg0zm3zf98Od7+37XRJqEDbe6qjnz3oWXzwg7ChZdt+S7MbCDeLpYqAOfcfc65yXib7scDxX7RQ6bPzMbg7fZagLcrZjCwHi+MwZs/x3Rg2iqBc8PmY1//e29Pt5bJDkxDuD8DX2qnTTcA4/B2bxyBtwuOkHF+grfrp9lRLQ117oBz7ofOuQl4xz7Ow1vucM4975w7Gy+I3vPbHT49ZwLfBb4KDPGn5+M2pqdNUZaFnXi/ysPXOc3fedTpCx11yPu2lpNOLxvxCII+ZtY35NXZ88MfBq4xsyn+0fIBZvYFMxvUzXY9Biwys+Fmlg38AO+AE3gHXoaZ2ZHNhc3sUjMb7v9irPU7N3agnpXA8WZ2iZllmNk/421aruhm+7tqKXCTmU0EMLMjzewr4YWcc43A74Fb/F9rE/COz7TlZbyVwyv+51XAv+Btgkf7rnYAR3d6KjqmM9/9m3j/jN8xsz5mNg34It5xkmazzewMM8vEO3D+pnOu0sxO8ZfPPv44PuXgshE+fQPw/qFrAMxsLt6v6Wb/ASw0s8n+8n6sv+INtxRY3NzPX44viMH3Ekl70xDuFuBM807lzvWHycbbf95sEN5B0FozGwr8e9g41gIX+/OmCO8HB/64Pm9mJ/hbD3vwVriNZhYws/P9H2z1eFuskZbDQXg/emqADDP7Ad4xhk6Ltiz4y//v8ObZIH++Xc/Bdc5a4CwzG+2vd25qp6oVwFFm9m9mluWPc4rfr9PLRjyCYCXeDG5+3dKZgZ1zpXhnUTyAt9/ufbz9e911G95ZQOuAt4G/+d1wzr2HFxQf+JtXI4FzgA1mtg/vQOHFzrlPO9D+XXi/UG7A20T+DnCec25nJ9v7Rzv0OoKnOjl8c3uewtuiWe5vgq8Hzo1SfAHebpbtePu6f9XO6F/G+6dqDoLX8H7lvBJ1CO/A2yL/e17YkWnoqM589865z4Dz8b6LncCDwOX+stDsN3grqN14ByS/7nc/Au8Hy0d4m/u7gLv9fo8AE/zpe9o59w7wE7x98juAE/DOBGlux+N4Z4f8Bu/g9NN4Wx7hfoZ3RtcLZrYX7+DglAjluvW9RBm+zWmIUL55n/ko4P/89v4v3pbJzX6xe/FODtjpT8tzYaO5Ge8X8Ed4x+p+E9LvKOAJvBB4F285/G+89dsNfj278XZRfitCE5/HO3a2CW/+fUrHdtdG0tay8C944fAB3v/Gb4BHAZxzL+KdLbkO72Bzm6HsnNsLnI33Y2U7sBn4vN+708uG+QcTREQkRR0OxwhERCSGFAQiIilOQSAikuIUBCIiKU5BICKS4hQEIiIpTkEgIpLiFAQiIilOQSAikuIUBCIiKU5BICKS4hQEIiIpTkEgIpLiFAQiIilOQSAikuIUBCIiKS6uQWBm55jZRjN738xujGfdIiISWdyeUOY/T3QT3uPVtgJrgK/5j70TEZFe0tkHyXfHqcD7zrkPAMxsOXABEDUIsrOzXX5+fnxaJyKSBMrKynY654Z3Zph4BkEuhz4QeivtPFA5Pz+f0tLSTlXinGP27NlceOGFXH311Z1u5O7du1m1ahVbt25tedXW1tLY2EhTUxONjY2dHqeISGcMGTKEJ598skvDmllFZ4eJZxBYhG6t9kuZ2VXAVQCjR4/ufCVmlJaW0tktib1793Lvvfdy9913s2fPHgCysrIYNWoUQ4YMIT09nbS0NNLT0zGLNCkiIj0j3j844xkEW4G8kM+jgKrwQs65h4CHAIqKirp0ACMQCLBjx44OlW1oaOD+++/n9ttvZ+fOnVx00UXccMMNjBs3jmHDhmmlLyJJL55nDa0BjjOzsWaWCVwM/CEWFXUmCO644w6uv/56TjrpJFavXs3vf/97Tj/9dLKzsxUCIpIS4rZF4JxrMLMFwPNAOvCoc25DLOoKBAKsXr263XJ79uzhnnvu4YILLuDpp5+ORVNERBJePHcN4ZxbCayMdT0d3SK4//77qa2t5Qc/+EGsmyQikrCS8sriQCDAvn372L9/f9Qye/fu5Z577uG8887j5JNPjmPrREQSS9IGAdDmVsGDDz7I7t27tTUgIikvJYNg37593H333Zx77rmccsop8WyaiEi7ysvL6devH4WFhRQWFnLNNdfEtL64HiOIl5ycHACqq6sj9l+6dCk7d+7k5ptvjmezREQ67JhjjmHt2rVxqSvltgj279/PkiVLOPvss/nc5z4X76aJSAopLy9n/PjxzJ8/n4kTJzJz5kzq6up6u1mtJGUQNG8RRAqCV199lerqaq6//vp4N0tEUtDmzZu59tpr2bBhA4MHD+bJJ59kyZIlLbt9Ql/XXXddy3BbtmzhpJNOYurUqbz66qsxbWNS7hrKyspi8ODBEYOgvLwcgIKCgji3SkRS0dixYyksLARg8uTJlJeXs2jRIoqLi6MOM2LECILBIMOGDaOsrIwLL7yQDRs2cMQRR8SkjUkZBBD9WoJgMEh6ejojRozohVaJSKrJyspqeZ+enk5dXR1LliyhpKSkVdmzzjqL++67j6ysrJbhJk+ezDHHHMOmTZsoKiqKSRtTMghGjRpFenp6L7RKRASKi4vb3CKoqalh6NChpKen88EHH7B582aOPvromLUnqYNg3bp1rboHg0HGjBnTCy0SEemYV155hR/84AdkZGSQnp7O0qVLGTp0aMzqS+ogiLZFcMYZZ/RCi0Qk1eTn57N+/fqWzwsXLuzQcHPmzGHOnDmxalYrSXnWEHhBUFtbS319fUu3xsZGtm7d2qXnHIiIJKukDgI49KKyDz/8kIaGBgWBiEiIpA+C0N1DwWAQ6NqTz0REkpWCQEQkxSkIRERSXMoFwZAhQxg0aFBvNUtEJOEkbRD069ePQYMGHRIEFRUV2hoQEQmTtEEAra8lCAaDCgIRkTAKAhGRFJfUQZCTk9MSBHv27KG2tlZBICISJqmDIHSLoLKyEkD3GRIRCdOjQWBmt5jZNjNb679mh/S7yczeN7ONZjarJ+uNJhAIsGvXLhoaGqioqAB06qiISLhY3HTup865u0M7mNkE4GJgIjASeMnMjnfONcag/hbNp5DW1NToGgIRkSjitWvoAmC5c67eObcFeB84NdaVhl5LEAwGycjI4Kijjop1tSIih5VYBMECM1tnZo+a2RC/Wy5QGVJmq98tpsKDQA+kERFprdNBYGYvmdn6CK8LgF8AxwCFwIfAT5oHizAqF2X8V5lZqZmV1tTUdLZ5hwgPAu0WEhFprdPHCJxzMzpSzsweBlb4H7cCeSG9RwFVUcb/EPAQQFFRUcSw6KjQIKioqGDq1KndGZ2ISFLq6bOGQp8IfxHQ/GiePwAXm1mWmY0FjgNW92TdkQwcOJB+/fpRVVXFtm3btEUgIhJBT581dJeZFeLt9ikHrgZwzm0ws98B7wANwLWxPmMIwMwIBAK89dZbNDY2KghERCLo0SBwzl3WRr/FwOKerK8jAoEAZWVlgE4dFRGJJKmvLAYvCPbt2wcoCEREIkmJIGimIBARaS1lgmDo0KEMHDiwl1sjIpJ4UiYItDUgIhKZgkBEJMUpCEREUlzSB0FOTg6gIBARiSbpg+DYY4/lW9/6FhdddFFvN0VEJCElfBDMmzePnJwcCgoK2i37yiuvcPLJJ5ORkcETTzwBQEZGBj//+c8ZN24chYWFFBYWcv7558e62SIihw1zrlv3dYspM6sBdgFNwFhgQzuDZALpQAD4GPgopN9JwFsxaKaISCIZ45wb3pkBEjoImplZPrDCOVfgfz4G+DkwHNgPzHfOvRdS/td++SdCuu1zzulCAhGRMAm/ayiKh4B/cc5NBhYCD3ZgmL7+cw7eMLMLY9o6EZHDSCyeWRxTZjYQOA143KzleTdZHRh0tHOuysyOBv5iZm875/4eq3aKiBwuDrsgwNuKqXXOFXZmIOdclf/3AzNbhXfMQEEgIinvsNs15JzbA2wxs68AmGdSW8OY2RAzy/LfZwOn4z0bQUQk5SX8wWIzewyYBmQDO4B/B/6C93zkEUAfYLlz7kdmdgrwFDAE+BTY7pybaGanAb/EO/soDbjXOfdIvKdFRCQRJXwQiIhIbPXIriEzO8fMNprZ+2Z2Y4T+Zmb3+f3XmdnJPVGviIh0X7eDwMzS8c7pPxeYAHzNzCaEFTsX74H1xwFX4e3WERGRBNATZw2dCrzvnPsAwMyWAxdw6MHYC4D/ct5+qDfMbLCZjXDOfdjWiLOzs11+fn6XGuWcI+T00pjS7jUR6WldXX+VlZXt7OyVxT0RBLlAZcjnrcCUDpTJBdoMgvz8fEpLSzvVGOcc/fv359vf/ja33357p4YFeOutt1i5ciXBYJDKykoqKyvZvXs39fX1La+Ghgat/EUkZgKBANu3b+/SsGZW0dlheiIIIsVW+FqyI2W8gmZX4e0+6tKto82MoUOHsmPHjk4N99prr7F48WKee+45AIYPH05eXh7HHHMMU6ZMoW/fvmRmZpKVlUVGRgZpaWmkpaVhZqSltb2HLV5bJiKSHAYMGBDX+noiCLYCeSGfRwFVXSgDgHPuIbxbSFBUVNSln92BQKDDQRAMBrnssst45ZVXyM7OZvHixXzzm99kyJAhXalaROSw0xNBsAY4zszGAtuAi4FLwsr8AVjgHz+YAnzc3vGB7uhMEBQXF1NWVsa9997L/Pnz6d+/f6yaJSKSkLodBM65BjNbADyPdwvoR51zG8zsGr//UmAlMBt4H+9uoXO7W29bAoEAGza0d8dqeOedd3j88ce58cYb+dd//ddYNklEJGH1yL2GnHMr8Vb2od2Whrx3wLU9UVdHNG8RtHfm0G233Ub//v25/vrr49U0EZGEc9jda6gjAoEAn332GR9//HHUMu+99x7Lly9nwYIFZGdnx7F1IiKJJWmDAGjzOMFtt91Gv379uOGGG+LVLBGRhJSSQbBp0yYee+wxrr32WoYP79R1FyIiPWb16tUtz1KfNGkSTz31VEu/adOmHfKs9erq6pi143B8HkG72guCxYsXk5WVpa0BEelVBQUFlJaWkpGRwYcffsikSZP44he/SEaGt2ouKSmhqKgo5u1IuS2CiooKSkpK+OY3v9lSTkSkO8rLyxk/fjzz589n4sSJzJw5k7q6unaH69+/f8tK/9NPP+21i0+Tcotg2LBhpKWlRQyC0tJSGhsbufTSS3uhZSKSrDZv3sxjjz3Gww8/zFe/+lWefPJJPvzwQ0pKSlqVPeuss7jvvvsAePPNN5k3bx4VFRUsW7asJRgA5s6dS3p6OnPmzGHRokUxC4qkDIL09HSGDx8eMQiCwSDg3cdIRKSnjB07lsLCQgAmT55MeXk5ixYtori4uM3hpkyZwoYNG3j33Xe54oorOPfcc+nbty8lJSXk5uayd+9e5syZw7Jly7j88stj0vak3DUE0a8uDgaDDBw4kMGDB8e/USKStLKyslrep6en09DQwJIlS1oO9oa+rrvuulbDjx8/ngEDBrB+/XoAcnNzARg0aBCXXHIJq1evjlnbk3KLAKIHQUVFBaNHj9aN4EQk5oqLi9vcItiyZQt5eXlkZGRQUVHBxo0byc/Pp6GhgdraWrKzszlw4AArVqxgxowZMWtn0gZBTk4OmzdvbtU9GAx26a6mIiI97bXXXuPOO++kT58+pKWl8eCDD5Kdnc0nn3zCrFmzOHDgAI2NjcyYMYP58+fHrB1JGwTRbjMRDAaZPHlyL7ZMRJJNfn5+yy4dgIULF3ZouMsuu4zLLrusVfcBAwZQVlbWY+1rT1IfI6irq2Pfvn0t3erq6qipqdEWgYhIiKQOAjj0WoLKSu8haWPGjOmVNomIJKKkD4LQy7KbTx3VFoGIyEFJHwShWwQKAhGR1lIuCMys5fxcERFJ4iBovqtoeBCMHDmSPn369FazREQSTrdOHzWzJcAXgc+AvwNznXO1EcqVA3uBRqDBORfz2+n16dOHYcOGtQoC7RYSETlUd7cIXgQKnHMnApuAm9oo+3nnXGE8QqBZ+NXFzVcVi4jIQd0KAufcC865Bv/jG8Co7jep54QGQVNTE5WVlQoCEZEwPXmMYB7wpyj9HPCCmZWZ2VU9WGebQoOgpqaG+vp6BYGISJh2jxGY2UvAURF6fd8594xf5vtAA9D6xtue051zVWaWA7xoZu85516JUt9VwFXQ/dM8Q4NAp46KiETWbhA459q85Z2ZXQGcB0x3zrko46jy/1ab2VPAqUDEIHDOPQQ8BFBUVBRxfB0VCATYu3cvdXV1CgIRkSi6tWvIzM4Bvguc75zbH6XMADMb1PwemAmsj1S2p4VeS6AgEBGJrLvHCB4ABuHt7llrZksBzGykma30ywSA18zs/4DVwLPOuee6WW+HhAfBwIEDGTJkSDyqFhE5bHTrOgLn3LFRulcBs/33HwCTulNPV4UHgR5IIyLSWtJeWQzew2ng0CAQEZFDKQhERFJcUgdB3759OfLIIykvL6e6ulpBICISQVIHAXjHCUpLSwGdMSQiEklKBEHzs0QVBCIiraVEEDQ2NgIKAhGRSFIiCAA9kEZEJIqUCYIRI0aQmZnZy60REUk8KRME2i0kIhKZgkBEJMWlTBCMGTOml1siIpKYkj4Img8Q5+fn925DREQSVMIHwbx588jJyaGgoKDdsvfccw8TJkzgxBNPZPr06VRUVJCXl8fTTz/N5s2bKSgooKCggN/+9rdxaLmIyOEh4YPgyiuv5LnnOnbX6pNOOonS0lLWrVvHl7/8Zb7zne8AkJGRwfr161m7di1vvvkmS5YsYc+ePbFstojIYcOiPFQsIZhZDVABZALHARv8XlnAaLzbaDf5ZT4NG7wfMAZ4D++ZCGnAh36/McAe4KMYNl9EpDeMcc4N78wACR0EzcwsH1jhnCvwP/8ZuMY5t9nMpgB3OOf+KWyYB4DtzrnbzGwm8O/A2UB/vAfk/Nw595N4ToeISCLq1oNpeoOZDQROAx4PechMVliZS4EiYCqAc+4FMzsFeB2oAf4KNMSrzSIiieywCwK8XTy1zrnCSD3NbAbwfWCqc66+ubtzbjGw2C/zG2Bz7JsqIpL4Ev5gcTjn3B5gi5l9BcA8k/z3JwG/BM53zlU3D2Nm6WY2zH9/InAi8ELcGy8ikoAS/hiBmT0GTAOygR14+/r/AvwCGAH0AZY7535kZi8BJ3DwoHDQOXe+mfUF/uZ324N3fGFt3CZCRCSBJXwQiIhIbB12u4ZERKRnxTUIzOwcM9toZu+b2Y3xrFtERCKL264hM0sHNuGdy78VWAN8zTn3TrRhsrOzne4RJCLScWVlZTs7e0FZPE8fPRV43zn3AYCZLQcuAKIGQX5+fsuD5zvKOceLL75Ibm4uEydO7E57aWxsZMeOHezevZv6+vqWV0NDA845mpqacM61vEREekJWVhYzZszo0rBmVtHZYeIZBLlAZcjnrcCU8EJmdhVwFXTtGQJmxkUXXcQ111zDT37SuQuH9+zZwy9+8QueffZZgsEg27Zto6FB152JSHwFAgG2b98et/riGQQWoVurn9HOuYeAhwCKioq69DM7EAiwY8eODpfftWsXP/vZz7j//vupra3l1FNP5cwzzyQvL4+8vDyys7Pp27cvmZmZZGVlkZGRQVpaGmaGmZGW5h1qCbnSWUSkyzIy4nutbzxr2wrkhXweBVTFoqKcnJwOB8GqVas477zz+OSTT/jSl77ETTfdRFFRUSyaJSKSkOIZBGuA48xsLLANuBi4JBYVBQIBtmzZ0m455xwLFy5k+PDhrF69mgkTJsSiOSIiCS1up4865xqABcDzwLvA75xzG9oeqms6umto5cqVlJWVsWjRIoWAiKSsuO6Ics6tBFbGup5AIMDOnTtpbGwkPT09Wlv44Q9/SH5+PpdffnmsmyQikrCS8sriQCBAU1MTu3btilrm+eefZ82aNXzve9+jT58+cWydiEhiSdogAKLuHmreGhg9ejRXXHFFPJsmInKIdevW8bnPfY6JEydywgkn8Omn3sMWp02bxrhx4ygsLKSwsJDq6up2xtR1h+PzCNoVGgQnnHBCq/4vvvgib7zxBkuXLiUzMzPezRMRAaChoYFLL72UZcuWMWnSJHbt2nXIHoqSkpK4nMWYclsEzVsDeXl5zJ07N95NE5EkVF5ezvjx45k/fz4TJ05k5syZ1NXVtTvcCy+8wIknnsikSZMAGDZsWNTjmrGUckHwt7/9jddff50bb7xRWwMi0mM2b97Mtddey4YNGxg8eDBPPvkkS5Ysadm1E/q67rrrANi0aRNmxqxZszj55JO56667Dhnn3LlzKSws5NZbb43pbWySctfQkUceSWZmZsQgePfddwGYPn16vJslIkls7NixFBYWAjB58mTKy8tZtGgRxcXFUYdpaGjgtddeY82aNfTv35/p06czefJkpk+fTklJCbm5uezdu5c5c+awbNmymJ3hmJRbBGYW9VqCYDAIQF5eXqt+IiJdlZWV1fI+PT2dhoaGdrcIRo0axdSpU8nOzqZ///7Mnj2bv/3Ne5hibm4uAIMGDeKSSy5h9erVMWt7Um4RQPSLyoLBYMuXLiISS8XFxW1uEcyaNYu77rqL/fv3k5mZycsvv8y3v/1tGhoaqK2tJTs7mwMHDrBixYou3420I5I6CKqqWt/KKBgMMmbMmF5okYjIoYYMGcL111/PKaecgpkxe/ZsvvCFL/DJJ58wa9YsDhw4QGNjIzNmzGD+/Pkxa0dCP7O4qKjIdfZ5BM2+8Y1v8Nxzz7Ft27ZDuhcUFHD88cfz+9//vieaKCKSUMyszDnXqXNOk/IYAXhbBNXV1TQ1NbV0c85RUVHRpecciIgkq6QOgoaGBj766KOWbh9//DH79u1TEIiIhEjqIIBDryVoPmNIQSAicpCCQEQkxSVtEOTk5AAKAhGR9iRtEETaIqioqCAzM7MlJEREJImDYOjQoaSnp7faIsjLy2t52LyIiCRxEKSlpbV6iH0wGNRuIRGRMEkbBND6NhMKAhGR1lImCA4cOEBVVZWCQEQkTMoEQVVVFU1NTbrPkIhImKQPgurqapxzOnVURCSKpA+C+vp69uzZoyAQEYki6YMAvGsJ9EAaEZHIUiYIKioq9EAaEZEIUiYIdOqoiEhkcQkCM1tiZu+Z2Toze8rMBsejXgWBiEj74rVF8CJQ4Jw7EdgE3BSPSrOzs0lLS2vZNaQgEBFpLS5B4Jx7wTnX4H98AxgVj3rT09PJzs5m48aNeiCNiEgUvXGMYB7wp3hVFggEWLNmDaBTR0VEIsnoqRGZ2UvAURF6fd8594xf5vtAA1DSxniuAq6CnllxBwIB3n77bXpqfCIiyabHgsA5N6Ot/mZ2BXAeMN0559oYz0PAQwBFRUVRy3VU6LMHFAQiIq31WBC0xczOAb4LTHXO7Y9Hnc2azxzKzMxseS8iIgfF6xjBA8Ag4EUzW2tmS+NUb8vKXw+kERGJLC5bBM65Y+NRTyTNQaDdQiIikSX9T2QFgYhI2xQEIiIpLumDID8/n379+jFp0qTeboqISEJK+CCYN28eOTk5FBQUdHiYJ554AjOjtLSUoUOH8te//pXbb7+dwsJCJk6cyNKlcTtWLSKS8BI+CK688kqee+65Dpffu3cv9913H1OmTGnpNn78eF5//XXWrl3Lm2++yZ133klVVVUsmisictixNq7t6nVmVgNUAJnAccAGv1cWMBrvrKcmv8ynfr88YA/eVc6VQPh1C+nABOA94EAMmy8i0hvGOOeGd2aAhA6CZmaWD6xwzhX4n/8MXOOc22xmU4A7nHP/ZGYnAYucc3PMbBWw0DlX6g+TBzwLHAsUO+d+3hvTIiKSaOJyHUFPMrOBwGnA42bW3DnLzNKAnwJXRhrOOVcJnGhmI4GnzewJ59yOODRZRCShHXZBgHdco9Y5Vxja0cyOBAqAVX5AHAX8wczOb94qAHDOVZnZBuBM4Im4tVpEJEEl/MHicM65PcAWM/sKgHkmOec+ds5lO+fynXP5eM89ON85V2pmo8ysn19+CHA6sLG3pkFEJJEkfBCY2WPAX4FxZrbVzL4BfB34hpn9H94B5AvaGc144E2//MvA3c65t2PZbhGRw8VhcbBYRERiJ65bBGZ2jpltNLP3zezGeNYtIiKRxW2LwMzS8R5cfzawFVgDfM05905cGiAiIhHF86yhU4H3nXMfAJjZcrx9+1GDIDs72+Xn58endSIiSaCsrGxnZy8oi2cQ5OJd6dtsKzAlSlnAu2FcaWlpW0Vacc5RXFzMmWeeyQUXtHcMObL9+/dTWVnZ8tq9ezf19fUtr4aGBpxzNDU14ZxreYmI9IRBgwaxePHiLg1rZhWdHSaeQWARurVae3b34fVmxiOPPEJ9fX2ngsA5x5///Gduv/12/ud//idquczMTDIyMkhLS8PMMLOWJ5+FXOAmItJlOTk5XQ6CrohnEGzFuw9Qs1FAqzu/9cTD6wOBADt2dOyiYeccf/zjH1m8eDGrV69m5MiR3HzzzRx//PHk5eWRl5dHdnY2ffv2pU+fPlrZi0jSiWcQrAGOM7OxwDbgYuCSWFTUmSB48MEHWbBgAWPHjuWXv/wlV1xxBVlZWbFolohIQopbEDjnGsxsAfA83h1AH3XObWhnsC4JBAKsW7eu3XJ1dXXcdtttTJ06lZdeeomMjMPxjhsiIt0T1zWfc24lsDLW9XR0i+Dhhx9m+/btLF++XCEgIikr4W8x0RWBQIDa2lrq6+ujlvn000/58Y9/zFlnncXUqVPj2DoRkcSSlD+Dmx9YX11dTV5eXsQyjzzyCFVVVSxbtiyeTRMRSThJu0UARN09VF9fz5133skZZ5zB5z//+Xg2TUSkxWeffcbcuXM54YQTmDRpEqtWrWrpN23aNMaNG0dhYSGFhYVUV1fHrB1JvUUQLQgeffRRtm7dyq9+9SudDioivebhhx8G4O2336a6uppzzz2XNWvWtFybVFJSQlFRUczbkZRbBDk5OUDkIKivr+eOO+7gtNNOY/r06fFumogkofLycsaPH8/8+fOZOHEiM2fOpK6urt3h3nnnnZb1UE5ODoMHD+703RR6QlIGQVtbBK+++iqVlZV897vf1daAiPSYzZs3c+2117JhwwYGDx7Mk08+yZIlS1p27YS+rrvuOgAmTZrEM888Q0NDA1u2bKGsrIzKyoN34pk7dy6FhYXceuutMb2NTVLuGurfvz8DBw6MGATl5eUAFBYWxrdRIpLUxo4d27JemTx5MuXl5SxatIji4uKow8ybN493332XoqIixowZw2mnndZyKntJSQm5ubns3buXOXPmsGzZMi6//PKYtD0pgwC8rYJIB1eCwSBpaWmMHDmyF1olIskq9I4E6enp1NXVsWTJEkpKSlqVPeuss7jvvvvIyMjgpz/9aUv30047jeOOOw6A3NxcwLsB3SWXXMLq1asVBJ0V7aKyYDBIbm6uLiATkZgrLi5uc4tg//79OOcYMGAAL774IhkZGUyYMIGGhgZqa2vJzs7mwIEDrFixghkzZsSsnUm7NgwEAmzatKlV92Aw2KW7moqI9LTq6mpmzZpFWloaubm5Ldc11dfXM2vWLA4cOEBjYyMzZsxg/vz5MWtHUgfBq6++2qp7MBjk1FNP7YUWiUiyys/PZ/369S2fFy5c2OHhNm7c2Kr7gAEDKCsr67H2tScpzxoCLwh27dpFQ0NDS7empiYqKyu1RSAiEiKpg8A5R01NTUu36upqPvvsMwWBiEiIpA4COPRagooK7wluCgIRkYNSKgiCwSCgIBARCaUgEBFJcSkXBIMGDeLII4/srWaJiCScpA2CgQMH0q9fv1ZBMGbMGN1jSEQkRNIGgZm1urpYF5OJiLSWtEEArW8zoSAQEWktZYJg//797Ny5U0EgIhImqYMgJyenJQia7/GtIBAROVRcgsDMvmJmG8ysycxi/9w1XyAQoKamhsbGRp06KiISRby2CNYDXwJeiVN9gBcETU1N7Nq1S1cVi4hEEZe7jzrn3gXiftpm6LUEeiCNiEhkCXeMwMyuMrNSMysNvWFcV4QHwciRI+nTp09PNFNEJGn02BaBmb0EHBWh1/edc890dDzOuYeAhwCKioq69bTm8CDQbiERkdZ6LAicc7F7jloXhQfBKaec0sstEhFJPAm3a6gnDR48mMzMTLZv364H0oiIRBGv00cvMrOtwOeAZ83s+TjVS05ODm+//TafffYZY8aMiUe1IiKHlXidNfQU8FQ86goXCARYs2YNoFNHRUQiSepdQ3Dw2cWgIBARiSQlgqCZgkBEpLWUCQI9kEZEJLKUCYLRo0frgTQiIhGkVBCIiEhrCgIRkRSnIBARSXFJHwRHH300U6dO5eyzz+7tpoiIJKS4XFDWm/r168eqVat6uxkiIgnLnOvWDT5jysxqgIreboeIyGFkjHNueGcGSOggEBGR2Ev6YwQiItI2BYGISIpTEIiIpDgFgYhIilMQiIikOAWBiEiKUxCIiKQ4BYGISIpTEIiIpLj/D+tdOz6/D8hTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "if __name__ == '__main__':\n",
    "    Program_3()\n",
    "    iflaga = 2\n",
    "    define_media(iflaga)\n",
    "    define_coefficients()\n",
    "    source = 2\n",
    "    ABC_order = 4;\n",
    "    define_Liao();\n",
    "    # Ez field plotting \n",
    "    for x in range(nt):\n",
    "        n = x+1\n",
    "        adv_Ez(n, source)\n",
    "        Liao_ABC()\n",
    "        adv_H(n)\n",
    "        #if n % 5 == 0:\n",
    "        #    my_surface_plot(Ez)\n",
    "        #    plt.title('3D plot of $E_z$ with PEC box')\n",
    "        #    plt.show()\n",
    "        if n == 5:\n",
    "            Ez_5 = copy.deepcopy(Ez[:, strip])\n",
    "        elif n == 35:\n",
    "            Ez_35 = copy.deepcopy(Ez[:, strip])\n",
    "        elif n == 65:\n",
    "            Ez_65 = copy.deepcopy(Ez[:, strip])\n",
    "        elif n == 95:\n",
    "            Ez_95 = copy.deepcopy(Ez[:, strip])\n",
    "\n",
    "    my_line_plot(Ez_5, Ez_35, Ez_65, Ez_95, iflaga, source)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8eed0875",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
